/*  sxmdcvdc.cpp  by K.Tsuru  */
// function ID = 504 BRADIX (reference)
/*****************************************************************************
SDecimal class
BRADIX ---> DRADIX
Provides the divide radix conversion.
Over 32,000 digits it is faster than ConvertToDec().But for 64,000 digits it
cannot be executed because of the out of memory,then not used in the library.
See
SLong SInteger::DConvToDec() const;
file name : simdcvdc.cpp.
******************************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
static uint SizeOfSDouble(){
    uint sz = sizeof(SDouble);
    sz += minArraySize*sizeof(fType);
    return sz+8u;
}
SDouble SDecimal::DConvToDec() const{
    if(Sign() == 0) return 0.0;
    uint N = Last()+1, rsz, z, s;
    int d, j, k, q;
    d = 16;
    rsz = min(N/d, SNMaxSize(REAL));
    if(!rsz) rsz++;     //divide number rsz >=1;   N = rsz*d - g; g >= 0
    uint unitsz = SizeOfSDouble();
    rsz =ceilpow2(rsz);
    while(rsz >= maxSizeOfMemoryBlock/unitsz +1){   // rsz is too large?
        d *= 2;     rsz /= 2;
    }
//evaluate again
    rsz = N/d;
    if(N % d) rsz++;    // rsz > 0

    RealSize C;
    C.SetEffFig(C.EffFigures()+C.Hidden(), C.TEMP_EXTEND);
    SDouble* r = new SDouble[rsz];

//Take into floating point mode for a moment.
    r[0].TempPointFree();
    SDouble R(1.0/(double)BRADIX), a;
        //1/BRADIX can be expressed exactly in double.
// step 1 :It converts into SDouble by d figures independently.
    for(k = 0; k < (int)rsz - 1 ; k++){ //pass if rsz == 1
        q = (k+1)*d - 1;
        r[k].SetInt((int)figure(q));    // Horner's method
        for(j = q - 1; j >= k*d; j--){
            a.SetInt((int)figure(j));
            r[k] = r[k]*R + a;  //faster than r[k] = DsDiv(r[k],BRADIX) + a;
        }
    }
//remaining upper part
    k = (int)rsz-1; r[k].SetInt((int)figure(N-1));
    for(j = (int)N-2; j >= k*d; j--){
        a.SetInt((int)figure(j));
        r[k] = r[k]*R + a;
    }
    r[0].PointModePop();    // r[0].CutDown(POP);
    if(rsz == 1) goto Ret;

//step 2 : It converts into SDouble binding two figures of r[].
//It is fast because the FFT multiplication is used.
    R.FixedPoint(1);
    z = ceilpow2(rsz);
    R = Dpow(R, d);         // = (1/BRADIX)^d

    while(z){
        uint i;
        for(i = 0; i < z ; i += 2){
            if(i+2u <= rsz) r[i/2] = r[i+1]*R + r[i]; //i <= rsz-2
            else if(i+1u == rsz) r[i/2] = r[i]; // i == rsz-1
            else r[i/2].SetZero(); // i >= rsz
        }
        z /= 2;
//free disused memory
        s = min(2*z, rsz);
        for(i = z; i < s; i++) r[i].figure.size(0, 0);
        if(z == 1) break;   //do not evaluate last R *= R;
        R *= R;
    }
    R.PointFree();
Ret:
//The last result is put in r[0].
    C.SetEffFig(0);
    R = r[0];
    delete[] r;
    if( R.Sign() ) R.SetSign(Sign());
    R.Reform(504);
    return R;
}